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I. INTRODUCTION 

Monte Carlo (MC) methods have become more pow- 
erful tools with the development of faster and more ac- 
cessible computers. Many different phenomena have been 
studied with MC methods . With standard MC (also 
sometimes called Metropolis or Markov Chain Monte 
Carlo) dynamics, one trial involves two parts; choosing a 
new state and deciding whether to accept or reject it. 

The standard dynamic MC procedure becomes very 
inefficient under some conditions, for example, at low 
temperature and in a strong external field. This is be- 
cause the rate of rejection becomes very high, so a huge 
number of trials is required to make a change in the state 
of the system. Various methods have been proposed to 
accelerate MC methods for studies of the statics of a sys- 
tem by modifying the underlying MC move P, Q . How- 
ever, when the underlying MC move is based on physical 
processes, such modifications of the MC method are not 
allowed since they would change the time development 
of the system. These kinetic MC methods are used in 
many physical situations, such as molecular beam epi- 
taxy m, as well as driven magnetic systems or models of 
membranes or biological evolution 

Accelerating MC methods without changing the un- 
derlying dynamics can be achieved using a different tech- 
nique called the rejection-free Monte Carlo (RFMC) 
method. The RFMC shares the original Markov chain 
with the standard MC, but it has rejection-less proce- 
dures. Therefore, a simulation of the RFMC is more 
efficient in the region where the standard MC is inefR- 
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cient due to many rejected trial states. The rejection- free 
scheme was first constructed for discrete spin systems , 
and has been applied for example to the kinetic Ising 
model in order to study dynamical critical behavior j^. 
For a review and history of the RFMC for discrete de- 
grees of freedom, see Ref. . A RFMC method has also 
been developed and applied to a model with continuous 
degrees of freedom 0. 

It is not trivial how to implement the rejection-free al- 
gorithm for each system. Therefore, it would be useful 
for a particular system to know how efficient the RFMC 
method is compared to the standard MC method without 
implementing a RFMC algorithm. The RFMC method 
has been applied to some spin systems. The standard 
MC method for particle systems can also become ineffi- 
cient in some conditions, for example, in the high den- 
sity or high dispersity. For example, it is important to 
study the nucleation and growth of defects such as the 
dislocations in a hard-disk system. While the dynam- 
ics of the defects are predicted by Kosterlitz-Thouless- 
Halperin-Nelson- Young theory ;8], there are few studies 
treating the nucleation because of the high-rejection rate. 
The hard-particle system with high dispersity is also of 
interest 9] . Such system can be a model of glassy mate- 
rials, and it is also difficult to study by the standard MC 
method. It is possible to use molecular dynamics (MD) 
simulations to study the time-dependent phenomena in- 
stead of MC. However, there are a number of difficulties 
also with using MD. For statics, both the MD and MC 
give comparable results (see for example 10] where ab 
initio MC and MD give comparable results and require 
comparable amounts of computer time). However, for 
dynamics neither the standard MC or the MD can go to 
long time scales, making studies, for example, of nucle- 
ation and growth computational unfeasible. Another dif- 
ficulty in MD simulations is that a particle-system has os- 
cillations of physical quantities because of the momentum 
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conservation. The oscihation cannot be removed by aver- 
aging independent samples, and prevents us from study- 
ing the dynamics of the order parameter Therefore, 
a rejection-free MC scheme for particle systems, as well 
as for spin systems, is also desirable. 

In this paper, we first give a brief review of the 
rejection-free scheme in Sec. ^ In Sec. IIIII we construct 
mean-field-type arguments which predict the efficiency of 
the RFMC compared to the standard MC. In Sec.|lV| we 
implement the RFMC method for the hard-disk system. 
Finally, we summarize our study and give discussions in 
Sec. El 



II. REJECTION-FREE SCHEME 



P(0|0) = A 




P(l|2) 



FIG. 1: A Markov cliain of Monte Carlo steps. 

A Monte Carlo method is an implementation of a 
Markov process on a computer, and hence is sometimes 
called a Markov Chain Monte Carlo. The Monte Carlo 
method calculates various physical quantities by updat- 
ing states of a system using random variables. These 
updating processes can be illustrated by a Markov chain 
(see the schematic in Fig. ^ . Let the current state be 
at 5*0. The states possible to move from Sq are denoted 
by Si (i = 1,2, ■ ■ ■ , Ng). Define Ei as the energy of the 
state Si. The new state Si {i — 0,1, - ■■ ,Ns) will be 
chosen with probability P{i\0). One way to ensure that 
the system will relax to the equilibrium state is to insist 
that the probability P{i\0) satisfies the detailed balance 
condition 0. 

One of the well-known ways to satisfy the detailed 
balance condition is to use a heat-bath transition prob- 
ability. In the heat-bath method, the probability P{i\0) 
is defined to be 



P(z|0) = 



exp {-PEi 



(1) 



When a system has a continuous degree of freedom, the 
summation of Eq. becomes an integration which is 
generally difficult to calculate analytically. 

Another popular way to satisfy the detailed balance 
condition is a Metropolis method. In this method, each 



step contains two parts; selecting a new state and accept- 
ing or rejecting the trial to move to the selected state. 
First, pick a state Si from all possible states to move to 
with uniform probability l/Ng. The probability P{i\Q) 
to move from 5*0 to Si is defined to be 1 when AEi < 0, 
otherwise it is exp {~f3AEi) with the energy difference 
AE., = El- Eq. Therefore, the probabihty P{i\0) in the 
Metropolis method is 



P(z|0) 



l/iV« 
exp {-(3AE,)/N, 



if AE, < 0, 
otherwise. 



(2) 



When a trial is rejected, the configuration of the system 
is not updated. The probability A = -P(0|0) to stay in 
the current state after the trial is given by 



A 



I-^P(^IO). 

i#0 



(3) 



For some parameters, e.g., under a strong external field 
and at an extremely low temperature, the value of A can 
be very nearly 1. In such cases, most of computational 
time is spent on calculating trials which will be rejected. 
This rejection rate drastically decreases the efficiency of 
the computation. For studies of the statics of a system, 
the MC trial move may be changed to increase the MC 
efficiency However, for studies of the dynamics 

a change in the MC move would change the underlying 
physics. 

In order to overcome this problem, a rejection-free 
Monte Carlo (RFMC) method is proposed. It is an ex- 
ample of an event driven algorithm UQ and has also 
been called a waiting time method [l^l • Each step of 
the RFMC method involves first computing the time to 
leave the current state (the waiting time i-^^rait)' ^^'^ then 
choosing a new state to move to with the appropriate 
probability. It does not contain the judgment to accept 
or reject a trial, and, therefore, it achieves rejection-less 
updates of the system in each algorithmic step. The wait- 
ing time i^a,it ^® ^ random variable. The probability p(t) 
to remain in the current state for t steps decays expo- 
nentially as, 



p{t) A* = exp (tin A), 



(4) 



with A defined in Eq. ©. Note that InA < since 
< A < 1. The time t to stay in the current state is 
determined to be, 



wait 



Inr 
ha 



1, 



(5) 



where r is a uniform random number on (0, 1) and \_x\ 
denotes the integer part of x. The rounding down is 
introduced to express the discrete time step in the MC 0, 
0. 

After the time of the system is advanced by ^-^^rait' ^ 
new state Si is chosen from the all states possible to move 
to, except for the current state, with the probability pro- 
portional to P{i\0) H H m. Since aU values of P{i\0) 
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are required to proceed one algorithm step in the RFMC, 
the computational cost of one step is higher than that of 
the normal MC. However, the waiting time t-^^ait' 
time which can be advanced in one algorithmic step, can 
become large, for example at low temperatures, and con- 
sequently the efficiency of the RFMC can become greater 
than that of the standard MC. 

It is worthwhile to stress that the RFMC method is 
mathematically equivalent to the standard MC method. 
Only the method of implementing the mathematics on a 
computer is different. Therefore, the dynamics are the 
same in both of the methods since they share the same 
Markov chain. This is in contrast to many other tech- 
niques to accelerate MC , which change the underlying 
relationship between single-trial MC time and the motion 
through the phase space of the system. 



III. EFFICIENCY OF RFMC 

In this section, we give arguments predicting the ef- 
ficiency of the RFMC method in spin and hard-particle 
systems. The efficiency of the RFMC method is inversely 
proportional to the rejection rate of the standard MC. 
Nevertheless, the RFMC method has the same dynamics 
as that of the standard MC method. Therefore, the effi- 
ciency of the RFMC is related directly to the inefficiency 
of the standard MC. 



A. Spin Systems 

Consider a general ferromagnetic spin system with MC 
dynamics at low temperature with a Hamiltonian, 



n 



SiS 



(6) 



with spins Si {\si\ = 1) and interaction energy J > 0. The 
sum is over the number of nearest-neighbor spins Ug ■ The 
expectation value of the waiting time is given by 



(''wait) 



1 



1-A' 



(7) 



which is the reciprocal of the acceptance probability . 
The rejection probability is A = Xi/N, with N the 
number of spins. The rejection probability of a MC move 
where spin i was the spin chosen to be changed is A^. 
At low enough temperature, the values of A^ are almost 
identical and any changes will usually involve an energy 
increase, since all of the spins are almost parallel. Ac- 
cordingly, the expectation value of Xi can be written as 



(A,) = 1 - (c-f^^^) 



(8) 



where /3 is the inverse temperature, AE is the energy 
difference between the current state and the chosen trial 



state and 



/sc 



denotes the average for all possible spin 



configurations. With Eqs. Q, and (|SJ, the expectation 
value of the waiting time is approximated by 



('wait) 



(9) 



The approximation involves replacing the expectation 
value of a function by the function of the expectation 
value, which is a mean-field or asymptotic type of ap- 
proximation. Equation implies that the waiting time, 
which is the efficiency of the RFMC method, is inversely 
proportional to the probability that the trial to flip a 
randomly chosen spin is accepted. Note that, the above 
argument depends only on the details of the spins, not 
on the lattice type or dimension of the system. In the 
following, we evaluate Eq. @ for three specific cases. 



1. Ising Model 

In the Ising model case, the energy difference AE is 
just Us J with a number of neighbor spins n^. The expec- 
tation value of the acceptance probability is (c~''^^)g(, = 



exp(— ris J/3), and therefore, the waiting time is 
('wait) ^ exp(nsJ/3), 



(10) 



which shows that the efficiency of the RFMC will in- 
crease exponentially as the temperature decreases. Simi- 
larly, for other systems with discrete degrees of freedom, 
such as the q-state Potts or clock models, i-^^^it increases 
exponentially with /3 (6|j. 



2. Classical XY Model 

When the spin has continuous degrees of freedom, the 
average in Eq. becomes an integration. For the clas- 
sical XY model, the expectation value of the acceptance 
probability becomes. 



1 



(11) 



since the energy increase AE = nsJ{l — cos6) with the 
angle of the spin 9. When (3^1, the integrand has a 
value only around 9 0. Therefore we make a saddle 
point approximation cos 9 ^ 1 — 0^/2, and change the 
upper limit of the integration to infinity. Then Eq. Illlj) 
reduces to the Gaussian integral. 



(e-^^")s 



1 



1 



With Eqs. © and we have. 



wait 



\/2imaJI3. 



(12) 



(13) 



Therefore, the efficiency of the RFMC method grows less 
rapidly with decreasing temperature in the XY model 
than it does for a discrete spin model. Nevertheless, at 
low enough temperature the RFMC will still outperform 
the standard dynamic MC. 
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3. Classical Heisenberg Model 

For the classical Heisenberg spin model, the energy dif- 
ference AE{9, (p) — ngj{\ — cos9), which is equivalent to 
the XY model. The expectation value of the acceptance 
probability is obtained from the integration, 



1 



2tt 



d6lsin6i-e-"=^''(^' 



2nsJP 



Therefore, 



(^wait) 2 



2ns J (3 



-2n,,JI3 



2nsJ(3, 



(14) 



(15) 



<C 1. This result, that the efficiency is pro- 
portional to /3, agrees with the past RFMC study 0- 
As the temperature is lowered, the efficiency of the 
RFMC for the classical Heisenberg spin model grows 
more rapidly than for the XY model, but not as rapidly 
as for a discrete spin model. 



B. Hard Particle Systems 



center of this chosen particle. The new position is cho- 
sen uniformly in the circle of radius CTs, called a step 
length, centered on the original position of the particle. 
This trial move is accepted when the new position has 
no overlap with any other particle, otherwise, the trial 
move is rejected. One MC step involves N trials, and 
the time evolution of this system can be considered to 
be Brownian-motion with a diffusion constant D oa (7^ 
for a low particle density. The probability Ai, which is 
the probability that the Monte Carlo trial of particle i is 
rejected, is 



A. -1 



A,. 



(16) 



with Ai the area particle i can move within ag with- 
out any overlap (see Fig. EJ. The rejection probabil- 
ity A is A = 1 — (A) /tt(t1, with the average of area 
{A) = J2i Ai/N. The mean distance between two neigh- 
boring particles is 2a and the radius of the particle is a. 
The radius of the area in which the particle can move is 
of order (a — cr) (see Fig.|3}. Therefore, 



A, 



(17) 



The density, p, of the system is inversely proportional 
to with the fixed radius cr, and p becomes the closest 
packing density pcp when a a. Therefore, we have 



Particles 




P 
Pep 



(18) 



From Eqs. ^ and ((TH|l . 

to be 

A,, 



the behavior of Ai is expected 



(19) 



with £ = (pcp — p)/ Pep, the result being valid for p near 
Pep- We thus obtain the expected value of the waiting 
time for the HD system to be, 



(^wait) ^ ^ £ 



(20) 



Similar arguments give {ty^^\^) for a d-dimensional hard 
particle system, 



FIG. 2: (Color online) A schematic drawing of the definition 
of Ai (shaded). The solid circles are particles and the small 
dashed circle has a radius cts. The shaded area is the area 
which is a continuous set of the points that the center of the 
chosen particle can move to. The ratio of Ai to the area of the 
trial circle ttcts^ gives the probability of accepting the move, 
1 — Ai, given that the center particle has been chosen as the 
one to move. 

Next, consider a hard-disk (HD) system with N par- 
ticles. A dynamic Monte Carlo procedure based on an 
underlying random-walk dynamics is: 1) choose one par- 
ticle randomly [from a uniform distribution over the in- 
dex of all N particles], 2) choose a new position for the 



(Vait) ^ [-7 



(21) 



Note that the behavior of (^wait) depends on the di- 
mension in the particle systems, while that of the spin 
systems does not. 



C. Simulations 

1. Waiting Time 

In order to confirm our predictions, Monte Carlo sim- 
ulations were performed on two- and three-dimensional 
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FIG. 3: (Color online) A schematic drawing of Ai. The 
radius of the particles and mean distance between centers of 
neighboring particles are denoted by a and 2a, respectively. 
The 'radius' of area Ai is on the order of (a — a), and hence 
Ai (X {a — a)^. 



Ising, XY, and Heisenberg spin systems on a square lat- 
tice (d = 2) and a simple cubic lattice (d = 3) . The linear 
system size simulated is L = 128, and periodic boundary 
conditions are used in all directions. After thermaliza- 
tion of 10'* MC steps per spin from the perfectly ordered 
state, the number of rejected trials Nr and the total tri- 
als Nt are counted over 10'^ MC steps per spin; there- 
fore, Nt = 128'^ ■ 10'^ with the dimension of the system 
d. Then the rejection probability A is approximated by 
A ^ Nr/Nt- Using this A, we estimated the value of 
(^wait)' "^^"^ simulation results are shown in Fig.^ The 
graphs show good agreement with our arguments. It is 
also worth noting that the prefactors we found are the 
same (within statistical errors) for the two- and three- 
dimensional systems. 

The behavior of the waiting time in the hard-particle 
systems is also confirmed. Monte Carlo simulations were 
performed on HD system with TV = 23288 and the hard- 
sphere (HS) system with N = 4000. After lO'* MC steps 
per particle, the value of {t^Q^[^} is estimated from 10* 
MC steps per particle. The results are shown in Fig. |S1 
in good agreement with our predictions. 
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2. Efficiency 

To further test our arguments, we implement the 
RFMC for the classical XY spin system. We discretize 
the spin state and use the 128-state clock model since we 
cannot calculate an acceptance probability analytically 
in this system. We confirmed that the behavior of the 
system with discretized spins is equivalent to the sys- 
tem with a continuum degree of freedom in the region 
where we simulated. At very low temperatures (lower 
than we simulated), where the number of states in the 
clock model approximation becomes important, we ex- 
pect that the RFMC efhciency crosses over to an expo- 
nential dependence as predicted in Eq. (^UJ). The sys- 
tem size is 128 x 128 and periodic boundary conditions 
are taken along the both directions. Measurements are 
started after 10^ MC steps. The CPU-time ratio of the 
standard MC to the RFMC methods to achieve 1000 ac- 



FIG. 4: (Color online) The waiting times t^a,it ^ '•'^ 
square lattices and simple-cubic lattices of (a) Ising, (b) XY, 
and (c) Heisenberg spin systems. Decimal logarithms are 
taken for the vertical axis of (a) and both axes of (b) and 
(c). Open circles are the calculated waiting time and solid 
lines are the predicted behavior with Ci = 0.48, C2 = 0.54, 
and C3 = 0.49. The number of neighboring spins Ws = 4 
for the two-dimensional and = 6 for the three-dimensional 
systems. They show excellent agreement with the predictions. 



cepted MC trials is shown in Fig. El The behavior of the 
efficiency of the RFMC is as predicted in Eq. (|13|l . 
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FIG. 5: (Color online) The waiting time t^^j^ of the hard- 
disk (HD) and hard-sphere (HS) systems at high densities. 
They are shown as functions of 1/e with e = (pep — p)/pcp. 
Decimal logarithms are taken for both axes. The solid lines 
are for visual reference with Ca = 0.07 and Cs — 0.04, respec- 
tively. This shows that the waiting time behaves as ~ e""^ 
with the dimension of the system d. 
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FIG. 6: (Color online) The ratio of the required computa- 
tional time to achieve 10'^ accepted MC trials of the standard 
MC to the RFMC method in the classical XY spin system. 
Decimal logarithms are taken along both axes. The solid line 
is for visual reference (C = 0.017). The efficiency of the 
RFMC behaves as predicted. 



IV. APPLICATION TO HARD-DISK SYSTEMS 

Consider a hard disk system with N particles aU with 
the radius a. Now that the high density expected ef- 
ficiency of the RFMC for hard particle algorithms has 
been obtained, an actual RFMC implementation for hard 
particles should be implemented. This section describes 
such an implementation. The standard MC method for 
the system involves choosing a particle, and trying to 
move the chosen particle within a circle with radius as 
centered on the original position of the chosen particle. 



To apply a RFMC method to the hard-disk system, de- 
fine Ai as the probability that a trial to move particle i 
is rejected (given that particle i was chosen as the par- 
ticle to attempt a move). Using the definition of A,;, we 
can construct the algorithm of the RFMC method for the 
hard-disk system as follows: 

1. Calculate the waiting time t-^g^n using Eq. (O with 

2. Advance the time of the system by i-^^rait- 

3. Choose a particle i with the probability propor- 
tional to 1 — Ai, which is the probability that (given 
that particle i was the particle chosen for an at- 
tempted move) the trial to move the particle i 
would be accepted. 

4. Choose the new position of the chosen particle i 
uniformly from all the points to which the particle 
i is allowed to move. 

The steps described above are the same as the RFMC 
for continuous spin systems 0, but the algorithms to 
calculate Ai, to choose a particle to move and to deter- 
mine a new position of the chosen particle are unique to 
the hard-disk system. In the following, we describe the 
details of the algorithms. 

A. Calculation of Ai 

The area Ai is the continuous set of positions in which 
the particle can be placed without any overlaps. Without 
neighboring particles, the shape of Ai would be a filled 
circle with a radius as- Let's call it a trial circle. In the 
general case, the shape of Ai is the remaining part of 
the trial circle after removing the overlap of 'shadows' of 
neighboring particles. The shape of the shadow is a circle 
with a radius 2a which is concentric to a neighboring 
particle. Let's call this a shadow circle. The area Ai, 
thus, consists of areas of arcs of a trial circle and that of 
shadow circles. 

To compute the value of Ai, we develop a method we 
call the survival point method. See Fig.|3(a). The cho- 
sen particle is shown as a solid circle, the trial circle is 
shown as a concentric dashed circle, and the area Ai is the 
shaded region. Each neighboring particle (filled circles) 
has a shadow circle which is concentric and has radius 2a. 
An enlargement of the area A, is shown in Fig. [7| (b). It 
is seen that in this example this area has five vertices 
which are intersection points of shadow particles, we call 
them survived vertex points. In Fig. [3(c), these survived 
vertex points are shown as small filled circles. Straight 
lines connect the center of the chosen particle and the 
intersection points. In this case the area Ai is divided 
into five portions. 

Each divided figure is the remaining part of an isosce- 
les triangle with the overlap of a shadow circle removed. 
It is easy to calculate this area. Thus, all we have to do 
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is to find all survived vertex points which form the area 
Ai. First, make a list of all intersection points of shadow 
circles and the trial circle. Next, remove points which 
are included in other shadow circles from the list, since 
these points cannot be vertices forming the area Ai. Af- 
ter this removal process, we have the vertices which form 
the area Ai (see Fig. [7|(c)). The calculation process of 
a partial figure which forms Ai is shown in Fig. {7\ (d). 
The vertices are denoted by Pi and P2, and the center 
of the shadow circle is denoted by S. The survived ver- 
tices Pi and P2 are on the shadow circle centered at S, 
so JP^ = 57^ = 2a. The area of OP1SP2 can be calcu- 
lated by summing the two triangles OP1P2 and SP1P2 
with Heron's formula. The area of the chord is Aa^O. 
Finally the portion of the area OP1P2 is calculated by 
subtracting the area of the chord SP1P2 from the area of 
the quadrilateral OP1SP2. The total area Ai is the sum 
of one such calculation for each survived vertex. 



B. Choosing a particle to move 

After calculation of t^Q_it and advancing the time of 
the system by it, we have to choose a particle i to move 
with a probability proportional to 1 — A^. With a direct 
implementation, i.e., with the integration scheme 
the order of the computation is 0{N), which is very time 
consuming. Other approaches are proposed like a three 
level search for spin systems The three-level search 
improves the efficiency of the search by determining co- 
ordinates of a spin to update one by one. However, it is 
difficult to apply this method for particle systems, since 
neighbors of particles are not fixed. Here we use a com- 
plete binary tree search for the choosing part of the al- 
gorithm. 

First, calculate the area Ai for each of the particles. 
Since an acceptance probability 1 — is proportional to 
Ai as shown in Eq. (|15|l . the particle should be chosen 
with the probability proportional to Ai . 

Next, construct a complete binary tree as follows, 

1 . Prepare a complete binary tree with enough height 
h, this height h should satisfy 2^"-^ < N < 2''~^ . 

2. Label each node with T^, which denotes the n}^ 
value at level k. The root node is labeled by T/*. 
A node labeled T^~^^ has branches leading to two 
nodes T^^_i and T^^. 

3. Associate every bottom node T.l with the value of 
area Ai. If the number of bottom nodes 2'*"^ is 
larger than N, the rest of the nodes are associated 
with zero, namely, = (z > N). 

4. Associate nodes at higher levels {k > 1) recursively 
with the sum of the values associated with its two 
children, namely, r,^i = r|„_i + T|„. 

A sample of a compete binary tree is shown in Fig. |51 
Each node has the value and the value of each node 




FIG. 7: (Color online) Computation of the value of Ai. (a) 
Filled gray circles represent neighboring particles with radius 
a and large circles are shadows of them (we call them 'shadow 
circles') with radius 2a. (b) An enlargement. The shaded area 
is the area into which the particle i is allowed to move. This 
figure is made by removing overlaps of shadow circles from 
the trial circle of radius as centered on the chosen particle, 
(c) The survived vertex points. The center of the trial circle 
is denoted by O. The solid circles represent survived vertex 
points, which form the area Ai. With them, we can calculate 
the value of Ai. The rectangle denotes a bounding rectangle. 
Each two adjacent survived vertex points and the center point 
O form a triangle. In this example there are five survived ver- 
tex points, and consequently five triangles to consider, (d) To 
calculate a portion of Ai, the area within each triangle formed 
by survived vertex points and O is calculated. The survived 
vertex points are the intersection points of the shadow circles 
or the trial circles, and here are denoted by Pi and P2. The 
center of the shadow particle is S. To find the shaded area 
a Monte Carlo procedure is performed in the shaded area of 
either (c) or (d), and only survived points generated in the 
shaded area are used as the new location for the new point 
O. 

at level fc -f 1 is the sum of the values of its two children 
nodes at level fc. The root node, which is in Fig. |S1 
has the sum of all Ai, that is, 

N 

Ti'^^^A,. (22) 

i 

Using this tree, we can choose a particle with the proba- 
bility proportional to Ai in the following way. 

1. fc ^ 1, i ^ 1. 

2. Prepare a random number r uniform on (0, T'/'). 

3 r z ^ 2i - 1 if r < T^-\ 
1 i <— 2i otherwise 
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A. k^i-1. 

5. if A; > 1 then go to|3 

Consequently, choosing the bottom node requires h — I 
random numbers. 

After the above processes, the final value of i indicates 
the index of the particle to move. The order of this search 
algorithm is O(logiV). When the position of particle i is 
moved, the value of Ai is also modified. Wc only have 
to update part of this tree for the chosen particle and 
its neighbors. The order of this update is also 0(log A^), 
which is much faster than 0{N) of the direct implemen- 
tation. Details to implement the complete binary tree 
search method are described in the appendix. 



Ti 



(t; 



Ti 



2n- 



1 + tL) 



Ti 



T{ Tj Ti Ti Ti Ti Tj Ti 



FIG. 8: A complete binary tree search. An example of A = 8 
[h — 4) is shown. The value of a node is the sum of the 



values of its two child nodes, namely, = TjT^ + ^2i-i- 
After construction of the tree, we use h — 1 random numbers 
to choose a bottom node. This bottom node is associated with 
the index of the particle that will be moved in this algorithmic 
step of the RFMC method. 



Carlo trial to find the new position of the particle to 
be moved became very inefficient. To improve this, it is 
effective to limit the trial area for the Monte Carlo by 
making the bounding area very close to the area Ai. We 
outline two different survived point methods, but have 
only implemented the first. 

For the first method, the one actually implemented 
in this paper see Fig. [71 (c). The solid rectangle is a 
bounding rectangle which includes the area Ai. It is 
easy to obtain the bounding rectangle with the survived 
vertex points. With the set of survived vertex points 
{(xi,yi)}, a diagonal line of the bounding rectangle is 
from (min {xi}, min {j/i}) to (max{a;i},max{2/i}). Then 
we can perform Monte Carlo trials for a new position 
within only this rectangle. The area of the rectangle is 
on the same order of Ai , so the probability of success to 
obtain the new position is drastically improved compared 
with the direct search over the trial circle. 

An alternative method is to first use a random number 
to decide which of the triangles formed with point O and 
two adjacent survived vertex points the survived point 
will fall into. This is done analytically since the areas 
of each triangle with removed shadow circle chords have 
been already calculated. Then the shortest side formed 
with point O and the two survived vertex points (say SP2 
in Fig. [71(d)) is lengthened to be equal to the longest side 
{SPi in Fig. [71(d)). The random trial point is then gen- 
erated within the section of the circle with a radius equal 
to the longest side {SPi in Fig. [3(d)). Then the point 
becomes the survived point used for the new location 
of point O if the trial point is within the shaded area. 
Otherwise, this procedure repeats in the same extended 
circular section until a survived point is found. 



D. Simulation 



C. Find a new position of the particle 



Calculation of Ai 



After choosing a particle, we have to choose the new 
position for it. It is difficult to choose a point uniformly 
from the points allowed to move into, since its shape is 
generally very comphcated (see Fig. [21). Therefore, we 
have chosen to choose the new position using a Monte 
Carlo rejection method. Namely we generate a random 
position uniformly over some bounding area that includes 
all of the area Ai. [Such as the dashed circle of radius CTs 
in Fig[7fa) or the rectangle in Fig [3c).] If this point is 
not in Ai it is rejected and another uniformly distributed 
random point over the bounding area is generated. The 
first point not to be rejected is the new position of the 
particle, since it is in the area Ai, and this point is used as 
the new center for the particle. Finally, the new Ai value 
for the chosen particle and all of its neighboring particles 
must be recalculated. This completes one algorithmic 
step of the RFMC method. 

The typical value of area Ai is very small compared 
to the trial circle at high density, and hence the Monte 



In order to test our method to calculate Ai described 
in Sec. IIV AI the values of Ai were also evaluated by a 
Monte Carlo sampling (Amc) with trial points uniformly 
drawn over the trial circle. The density of the system p is 
defined to be p = ANa'^/L'^ with the number of particles 
A, the radius of the particles a and the linear system 
size L, respectively. Throughout this study, the number 
of particles N is set to be 23288 and periodic bound- 
ary conditions are taken for both axes. The number of 
the generated configurations were 3000, and 10^ MC trial 
points are taken for each of the configurations to evaluate 
its area. The density of the system is fixed at p = 0.9. 
The result is shown in Fig. [HI The area Ai is normalized 
by the area of the trial circles (see Fig.l^J. The difference 
between the MC and our survived point method is less 
than 0.01% for all areas, which is the same order as the 
statistical error of our MC method. The standard devia- 
tion of the MC area calculation is determined by dividing 
the data into 10 groups, each including 10^ samples. This 



result shows that the value of Ai is properly calculated 
by our method. 




0.98 



0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

AMc/(7ro-r) 

FIG. 9: Comparison of calculated Ai between our survived 
points method and that calculated by a more straightforward 
MC method. Units on both axes are dimensionless. The cir- 
cles denote the ratio of j4-p^p|y[Q to the A-^q. The number 
of configurations is 3000 at p = 0.9 and 10*' independent sam- 
ples are averaged for each configuration. The solid lines are 
the standard deviation obtained from the jack-knife procedure 
described in the text. 



2. Time evolution 

The dynamics of the standard MC and of the RFMC 
must be the same. To test this in our case, we observed 
the time evolution of the six- fold bond-orientational order 
parameter 06 8]. The parameter 06 is defined to be 



06 = (exp(6«6')) , 



(23) 



with the bond angle 6 which has a definition described in 
Fig. ^1 The average is taken for all pairs of neighboring 
particles. The parameter 06 becomes 1 when all parti- 
cles are located on the points of a hexagonal grid, and it 
becomes when the particle location is completely disor- 
dered. Therefore 06 describes how close the system is to 
the perfect hexagonal packing. The neighbors in an off- 
lattice model are strictly defined with the Voronoi con- 
struction [T^ . which is a very time-consuming method. 
In this paper, two particles separated by a distance less 
than 2.6a are defined as neighbors. We confirmed that 
the value of 06 is approximately the same value as the 
value obtained with the Voronoi construction. At the 
beginning of the simulation, the particles are set up in a 
perfect hexagonal order, namely, 06(t — 0,p) = 1. The 
order parameter 06 starts to relax to the value of the 
equilibrium state. With this nonequilibrium relaxation 
(NER) behavior of order parameters, critical points and 
critical exponents of vario us p hase transitions can be de- 
termined accurately |3 l2Ct l2ll |2^ . This method is 
called a NER method. Watanabe et al. [H ^ stud- 
ied two-dimensional melting based on the NER method 




Area to search 



Definition of the angle 



FIG. 10: The definition of the neighboring particles and bond 
angle 9. Two particles separated by a distance less than R 
are defined as neighbors. Here, R is set to be 2.6a with the 
radius a of particles. The bond angle 9 is defined to be an 
angle between the bond connecting neighboring particles with 
respect to an arbitrary, but fixed global axis. 



for the Kosterlitz-Thouless transition |25| by observing 
the relaxation behavior of 06. Therefore, the following 
time evolutions of 06 contains information about the two- 
dimensional melting transition. 

Time evolutions of 06 are shown in Fig. 1111 Solid 
lines are results of the standard MC simulation and sym- 
bols (circles, triangles and squares) are the results of the 
RFMC. Fig.^Jshows that both behaviors are equivalent 
for the two methods. This is essentially a check of the 
program implementation, since the physical dynamic is 
the same for both the MC and the RFMC methods. 
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FIG. 11: Relaxation behaviors of the bond-orientational or- 
der. Solid lines are the results of the standard MC and sym- 
bols (circles, triangles and squares) are the results of the 
RFMC. The data intervals between accepted updates for the 
RFMC algorithm becomes longer at high density, while the 
data keeps the behavior of the standard MC algorithm. 



3. Efficiency 

The computational times required to achieve 1000 ac- 
cepted MC steps are shown in Fig.^J^a). Configurations 
are started from the perfect hexagonal configuration and 
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FIG. 12: (a) The required computational time to achieve 1000 
acceptances of the Monte Carlo moves with the standard MC 
(open circles) and the RFMC (solid circles), (b) CPU-time 
ratio vs. 1/e with e = (pep ~ p)/pcp. Decimal logarithms are 
taken for both axes. The solid line is drawn for the visual 
reference (C = 0.9 x 10~^). 



both measurements are started after 10^ MC steps. All 
simulations are performed on an Intel Xeon 2.4 GHz com- 
puter. While the computational time of the RFMC (open 
circles) is almost constant, a longer computational time 
is required for the standard Monte Carlo (solid circles) at 
higher density. It shows that the RFMC is more efficient 
at high densities, in spite of the additional bookkeeping 
involved in the RFMC method (so one RFMC algorith- 
mic step takes much longer than one standard MC step). 

The CPU-time ratio of the standard MC to the RFMC 
is shown in Fig. ll2r b). The data are shown as a function 
of 1/e, where e = (pcp — p)/pcp, the closest density is 
Pep and the density of the system is p. The CPU-time 
ratio, which is the efficiency of the RFMC compared to 
that of the standard MC, diverges as £~^. This confirms 
the predicted behavior of Eq. (|21|) . 



V. SUMMARY AND DISCUSSION 

We predicted the behavior of the average waiting time 
(*wait) to be 



(^wait) 



exp (const. /3) 



Ising and other discrete 
spin systems 
^ (classical XY) 
/3 (classical Heisenberg) 

(hard-particle) 

(24) 

for the efficiency of the RFMC method. These have 
been confirmed by our MC simulations. It is interesting 
that the behavior of the HD system in the high density 
regime is different from that of the XY spin system at 
low temperature, while the phase transition of both are 
of a Kosterlitz-Thouless-type 8, 26, 27] (see for a review 
Ref. 13). Our arguments for (t-^^ait) ^^^^ general, 
and consequently should be able to give the RFMC effi- 
ciency for other models, for example, a discrete stochastic 
model 2%. 

We implemented the rejection-free Monte Carlo algo- 
rithm for the hard-disk system. This method conserves 
the property of the dynamic behavior of the original 
Monte Carlo method. In other words, the time scales 
will not depend on the density, but are rather set by some 
Brownian-motion type of dynamic for all densities. An 
estimate of the time scales between the MC and physical 
time can thus be obtained by setting the mean-free path 
of an isolated particle to be proportional to the value Os ■ 
Note that strictly this is only true in the limit ~> 0, 
but it should be a reasonable approximation for a small 
finite Us- 

We also find that for a fixed value of CTs, the RFMC 
method is more efficient at high density. Therefore, 
the RFMC method should be useful for studies of two- 
dimensional solids or studies of high-density glass mate- 
rials, while the efficiency of the RFMC is less than that 
of the standard method at the critical point. It may 
also be possible to make the algorithm even more effi- 
cient by further optimization techniques, e.^., by utilizing 
the ideas of absorbing Markov chains (for the MCAMC 
method for discrete state spaces see Q and references 
therein) . Increased algorithmic efficiencies for the Monte 
Carlo dynamics of hard disks could be useful to further 
test physical phenomena using hard disk systems, such 
as for example the relationship between fluctuations and 
dissipation of work in a Joule experiment [sO] . 

The RFMC method gives the same dynamics as the 
standard MC method, and consequently, the RFMC 
method allows one in certain regimes to efficiently study 
the dynamical behavior of systems with a given physi- 
cal MC dynamic. The dynamic for a MC has been de- 
rived from underlying physical properties for some sys- 
tems 1^ Is^, • It has been shown that using differ- 
ent MC dynamics can have enormous consequences on 
physic al prop erties such as on low-temperature nucle- 
ation |34l l35| . Consequently, this equivalence between 
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the two MC methods is essential. 
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i ^ N 
While i ^ 

a{i) <- a{2i) + a{2i + 1) 

I <— i — 1 
next i 

Using this array, we can pick particle i with the proba- 
bility proportional to Ai as. 



While i < N 

Prepare a uniform random number r on (0, a{i)) 



i <— 2i if r < a{2i) 

i -i— 2i + l otherwise 



Appendix 



a(l) a(2) a(3) ffl(4) a(5) a(6) a(7) a.{N) a(2N - 1) 



next i 
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FIG. 13: Implementation of the complete binary tree search 
with an array. The required size of the array to implement 
the tree with height h is 2^ ~ \. The height h should satisfy 

2^-^ < N < 2^'^. When the number of particles iV is 2^'^, 
which is the maximum number of particles that the tree with 
height h can treat, the size of the array is 2N — 1. 



The complete binary tree search can be implemented 
with an one-dimensional array. To make it simple, let 
the number of particles N be 2^*"^. The tree with height 
h requires an array a(i) with size 2N — 1. First, associate 
each bottom node with a corresponding value as 



a{N + Ai (i = l,2,---,7V), 



(25) 



which corresponds to <— Ai. Next, associate parent 
nodes recursively as 



N+1. 



After the above procedure, we obtain the index i of the 
chosen particle. When the value of Ai is changed, the tree 
should be updated. The update process is as follows, 



a{N + 

i + N 



While i ^ 1 

a{i) ^ a{2i) + a{2i + 1) 
i ^ [i/2\ 
next i. 

Note that, when the chosen particle is moved, the ac- 
ceptance probabilities of the neighboring particles of the 
moved particle are also changed. Therefore, we have to 
perform the above process for all neighboring particles. 



[1] K. Binder, Monte Carlo Methods in Statistical Physics, 
ed. K. Binder, (Springer, Berlin, 1979). 

[2] D. Landau and K. Binder, A Guide to Monte Carlo Simu- 
lations in Statistical Physics, 2nd Ed. (Cambridge, 2005). 

[3] J. Hausmann and P. Rujan, Phys. Rev. Lett. 79, 
3339 (1997); P.B. Sunil Kumar, G. Gompper, and R. 
Lipowsky, Phys. Rev. Lett. 86, 3911 (2001); D. Chowd- 
hury, D. StaufFer, and A. Kunwar, Phys. Rev. Lett. 90, 
068101 (2003). 

[4] A. B. Bortz, M. H. Kalos and J. J. Lebowitz, J. Comput. 

Phys 17, 10 (1975). 
[5] S. Miyashita and H. Takano, Prog. Theor. Phys. 63 797 



(1980). 

[6] M. A. Novotny in Annual Reviews of Computational 
Physics IX, editor D. Stauffer (World Scientific, Singa- 
pore, 2001), p. 153. 

[7] J. D. Mufioz, M. A. Novotny and S. J. Mitchell, Phys. 
Rev. E 67, 026101 (2003). 

[8] B. I. Halpcrin and D. R. Nelson, Phys. Rev. Lett. 41, 
121 (1978);A. P. Young, Phys. Rev. B 19, 1855 (1979). 

[9] H. Watanabe, S. Yukawa, and N. Ito, Phys. Rev. E, 71, 
016702 (2005). 

[10] S. Wang, S.J. Mitchell, and P.A. Rikvold, Comp. Mater. 
Sci., 29, 145 (2004). 



12 



[11] H. Watanabe, S. Yukawa, and N. Ito, submitted to Phys. 
Rev. E. 

[12] J. Liu and E. Luijtcn, Phys. Rev. Lett. 92 035504 (2004). 
[13] J. Ball and P. Sibani, Comp. Phys. Commun. 141 260 
(2001). 

[14] M. A. Novotny, Phys. Rev. Lett. 74, 1 (1995); erratum 
75 1424 (1995). 

[15] T. Koma, J. Stat. Phys. 71, 269 (1993). 

[16] The order of the search itself is 0(lnA'') in the direct 
implementation with the integration scheme, since we can 
use a binary search for an ordered list. However, when one 
value is modified, the order of the computation to update 
the list of values is 0{N). 

[17] J. L. Blue, I. Beichl, and F. Sullivan, Phys. Rev. E 51, 
R867 (1995). 

[18] A. Jaster, Phys. Rev. E 59, 2594 (1999). 

[19] N. Ito, Physica A 192 (1993) 604. 

[20] N. Ito, Physica A 196 (1993) 591. 

[21] N. Ito, T. Matsuhisa and H. Kitatarii, J. Phys. See. Jpn. 
67 (1998) 1188. 

[22] N. Ito, K. Hukushima, K. Ogawa and Y. Ozeki, J. Phys. 

See. Jpn. 69 1931 (2000). 
[23] H. Watanabe, S. Yukawa, Y. Ozeki and N. Ito, Phys. 

Rev. E 66, 041110 (2002). 



[24] H. Watanabe, S. Yukawa, Y. Ozeki, and N. Ito, Phys. 

Rev. E, 69, 045103(R) (2004). 
[25] Y. Ozeki, K. Ogawa and N. Ito, Phys. Rev. E, 67, 026702 

(2003). 

[26] H. Watanabe, S. Yukawa, Y. Ozeki, and N. Ito, Phys. 

Rev. E 66, 041110 (2002); 
[27] Y. Ozeki, K. Ogawa, and N. Ito, Phys. Rev. E 67, 026702 

(2003); Y. Ozeki and N. Ito, Phys. Rev. B 68, 054414 

(2003). 

[28] K.J. Strandburg, Rev. Mod. Phys. 60, 161 (1988). 
[29] B.T. Werner and D.T. Gillespie, Phys. Rev. Lett. 71, 
3230 (1993). 

[30] B. Cleuren, C. van den Broeck, and R. Kawai, Phys. Rev. 

Lett. 96, 050601 (2006). 
[31] P.A. Martin, J. Stat. Phys. 4, 194 (1977). 
[32] K. Park, M.A. Novotny, and P.A. Rikvold, Phys. Rev. E 

66, 056101 (2002). 
[33] X. Z. Cheng, M.B.A. Jalil, H.K. Lee and Y. Okabe, Phys. 

Rev. Lett. 96, 067208 (2006). 
[34] K. Park, P. A. Rikvold, G. M. Buendi'a, and M. A. 

Novotny Phys. Rev. Lett. 92, 015701 (2004). 
[35] G. M. Buendi'a, P. A. Rikvold, K. Park, and M. A. 

Novotny, J. Chem. Phys. 121, 4193 (2004). 



